9. Diseño de sistemas y filtros IIR¶

Un filtro FIR de buena calidad puede requerir una gran cantidad de coeficientes

Es posible implementar filtros más eficientes usando recursividad. Esta es la base de los filtros de respuesta al impulso infinita o IIR que veremos en esta lección

9.1. Definición de un sistema IIR¶

Generalizando el sistema FIR para incluir versiones pasadas de la salida y asumiendo \(a[0] = 1\) llegamos a

\[\begin{split} \begin{align} y[n] &= b[0] x[n] + b[1] x[n-1] + b[2] x[n-2] + \ldots + b[L] x[n-L] \nonumber \\ & - a[1] y[n-1] - a[2] y[n-2] - \ldots - a[M] y[n-M] \nonumber \\ &= \sum_{l=0}^{L} b[l] x[n-l] - \sum_{m=1}^{M} a[m] y[n-m] \nonumber \\ \sum_{m=0}^{M} a[m] y[n-m] &= \sum_{l=0}^{L} b[l] x[n-l] \nonumber \\ (a * y)[n] &= (b * x)[n], \nonumber \end{align} \end{split}\]

es decir dos convoluciones discretas que definen una ecuación de diferencias

Este tipo de sistema se conoce como

sistema infinite impulse response (IIR)
sistema auto-regresive moving average (ARMA)
- autoregresivo de orden M: incluye valores pasados de la salida
- media movil de orden L+1: pondera el valor presente y pasados de la entrada

Podemos ver el sistema IIR como una generalización del sistema FIR. El caso particular del sistema FIR se recupera si

\(a[m] = 0\) para \(m=[1, \ldots, M]\)

9.1.1. Respuesta en frecuencia del sistema IIR¶

Aplicando la transformada de Fourier convertimos las convoluciones en multiplicaciones y encontramos la respuesta en frecuencia como

\[\begin{split} \begin{align} \text{DFT}_N[(a * y)[n]] &= \text{DFT}_N[(b * x)[n]] \nonumber \\ A[k] Y[k] &= B[k] X[k] \nonumber \\ H[k] = \frac{Y[k]}{X[k]} &= \frac{B[k]}{A[k]} = \frac{ \sum_{l=0}^L b[l]e^{-j \frac{2\pi}{N} nl} }{ \sum_{m=0}^M a[m]e^{-j \frac{2\pi}{N} mk}} \nonumber \end{align} \end{split}\]

que existe siempre que \(A[k] \neq 0\).

La respuesta en frecuencia también suele expresarse como

\[ H[k] = K \frac{ \prod_{l=1}^L (e^{j \frac{2\pi}{N} k}- \beta[l]) }{ \prod_{m=1}^M (e^{j \frac{2\pi}{N} k}- \alpha[m])} \]

donde

\(K\) se llama ganancia
las raices del polinomio del numerador \(\alpha\) se llaman conjuntamente ceros
las raices del polinomio del denominador \(\beta\) se llaman conjuntamente polos

9.1.2. Ejemplo de respuesta al impulso de un sistema IIR¶

Consideremos el siguiente sistema IIR

\[\begin{split} \begin{align} y[n] &= (1-\gamma) x[n] + \gamma y[n-1] \nonumber \\ y[n] - \gamma y[n-1] &= (1-\gamma) x[n] \nonumber \end{align} \end{split}\]

Los coeficientes del sistema son

\(a[0] = 1\), \(a[1] = -\gamma\) y \(b[0] = (1-\gamma)\)

Es decir que es AR de orden 1 y MA de orden 1

¿Cúal es su respuesta al impulso? Asumiendo \(y[n]=0, n<0\), tenemos que

\[\begin{split} \begin{matrix} n & \delta[n] & y[n] \\ -2 & 0 & 0 \\ -1 & 0 & 0 \\ 0 & 1 & (1-\gamma) \\ 1 & 0 & \gamma(1-\gamma) \\ 2 & 0 & \gamma^2(1-\gamma) \\ 3 & 0 & \gamma^3(1-\gamma) \\ 4 & 0 & \gamma^4(1-\gamma) \\ \end{matrix} \end{split}\]

¿Cómo cambia la respuesta al impulso con distintos valores de \(\gamma\)? ¿Qué pasa si \(\gamma \geq 1\)?

Respondamos estas preguntas visualizando la respuesta al impulso de este sistema con la función scipy.signal.dimpulse

# Valores de gamma que probaremos:
gamma = [-1.5, -1, -0.5, 0.5, 1., 1.5]

p = []
for g in gamma:
    t, y = scipy.signal.dimpulse(([1-g, 0], [1,-g], 1), x0=0, n=30)
    p.append(hv.Curve((t, y[0][:, 0]), label=f"gamma={g}"))
    
hv.Layout(p).cols(3).opts(hv.opts.Curve(width=250, height=200, axiswise=True))

/home/phuijse/.conda/envs/info320/lib/python3.9/site-packages/scipy/signal/filter_design.py:1630: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless
  warnings.warn("Badly conditioned filter coefficients (numerator): the "

De las figuras podemos ver que:

Para \(\gamma < 0\) (primera fila) los coeficientes del sistema son alternantes en signo
Para \(|\gamma| < 1\) los coeficientes del sistema tienden a cero
Para \(|\gamma| > 1\) los coeficientes del sistema divergen y tienen a infinito

Advertencia

A diferencia de un sistema FIR, el sistema IIR puede tener configuraciones inestables en que los coeficientes crecen o decrecen infinitamente

Por otro lado consideremos el sistema anterior y asumamos que \(|\gamma|<1\), desenrollando tenemos que

\[\begin{split} \begin{align} y[0] &= (1-\gamma) x[0] \nonumber \\ y[1] &= (1-\gamma) (x[1] + \gamma x[0]) \nonumber \\ y[2] &= (1-\gamma) (x[2] + \gamma x[1] + \gamma^2 x[0]) \nonumber \\ y[3] &= (1-\gamma) (x[3] + \gamma x[2] + \gamma^2 x[1] + \gamma^3 x[0]) \nonumber \\ y[4] &= (1-\gamma) (x[4] + \gamma x[3] + \gamma^2 x[2] + \gamma^3 x[1] + \gamma^4 x[0]) \nonumber \\ y[5] &= \ldots \nonumber \end{align} \end{split}\]

Nota

Con un sistema IIR de pocos coeficientes podemos representar un sistema FIR considerablemente más grande

En el ejemplo anterior, si escogemos \(\gamma\) tal que \(\gamma^{20 }\approx 0\) entonces aproximamos un sistema FIR de orden 20 con tan sólo 3 coeficientes

9.1.3. Ejemplo de respuesta en frecuencia de un sistema IIR¶

Para el sistema del ejemplo anterior su respuesta en frecuencia es

\[\begin{split} \begin{align} Y[k] &= (1-\gamma) X[k] + \gamma Y[k] e^{-j \frac{2\pi}{N} k} \nonumber \\ H[k] = \frac{Y[k]}{X[k]} &= \frac{1-\gamma}{1 - \gamma e^{-j \frac{2\pi}{N} k} } \nonumber \end{align} \end{split}\]

que en notación de polos y ceros se escribe como

\[ H[k] = (1-\gamma)\frac{e^{j \frac{2\pi}{N} k} - 0}{e^{j \frac{2\pi}{N} k} - \gamma } \]

es decir que tiene un cero en \(0\), un polo en \(\gamma\) y una ganancia de \((1-\gamma)\)

Para entender mejor este sistema estudiemos la magnitud de \(|H[k]|\) para \(\gamma < 1\)

\[\begin{split} \begin{align} | H[k]| &= \frac{|1-\gamma|}{|1 - \gamma e^{-j \frac{2\pi}{N} k}|} \nonumber \\ &= \frac{1-\gamma}{\sqrt{1 - 2\gamma \cos(\frac{2\pi}{N} k) + \gamma^2}} \nonumber \end{align} \end{split}\]

¿Cómo se ve \(|H[k]|\)? ¿Qué función cumple este sistema?

k = np.arange(-24, 25)/50
Hk = lambda gamma, k : (1-gamma)/np.sqrt(1 - 2*gamma*np.cos(2.0*np.pi*k) + gamma**2)

p = []
for gamma in [0.25, 0.5, 0.75]:
    p.append(hv.Curve((k, Hk(gamma, k)), 'Frecuencia', 'Respuesta', label=f'gamma={gamma}'))
    
hv.Overlay(p)

Nota

Este sistema atenua las frecuencias altas, es decir que actua como un filtro pasa bajos

9.2. Diseño de filtros IIR simples¶

Los filtros IIR más simples son los de un polo y un cero, es decir filtros de primer orden

\[ H[k] = \frac{b[0] + b[1] e^{-j \frac{2\pi}{N} k}}{1 + a[1] e^{-j \frac{2\pi}{N} k}} = K\frac{e^{j \frac{2\pi}{N} k} - \beta}{e^{j \frac{2\pi}{N} k} - \alpha } \]

donde podemos reconocer

\(b[0]=K\)
\(\beta = - b[1] \cdot K\)
\(\alpha=-a[1]\)

Definimos la frecuencia de corte \(f_c\) como aquella frecuencia en la que el filtro alcanza una atenuación de 0.7 (-3 dB). Haciendo la equivalencia con el ejemplo anterior tenemos que \(\gamma = e^{-2\pi f_c}\)

9.2.1. Receta para un filtro pasa bajo IIR con frecuencia de corte \(f_c\)¶

Asignamos

\(b[0] = 1 - e^{-2\pi f_c}\)
\(b[1] = 0\)
\(a[1] = -e^{-2\pi f_c}\)

Lo que resulta en la siguiente respuesta en frecuencia

\[ H[k] = \frac{1-e^{-2\pi f_c}}{1 - e^{-2\pi f_c} e^{-j \frac{2\pi}{N} k}} = (1-e^{-2\pi f_c}) \frac{(e^{j \frac{2\pi}{N} k}- 0)}{(e^{j \frac{2\pi}{N} k} - e^{-2\pi f_c} )} \]

Es decir un cero en \(0\), un polo en \(e^{-2\pi f_c}\) y ganancia \(1-e^{-2\pi f_c}\)

9.2.2. Receta para un filtro pasa alto IIR con frecuencia de corte \(f_c\)¶

Asignamos

\(b[0] = (1 + e^{-2\pi f_c})/2\)
\(b[1] = -(1 + e^{-2\pi f_c})/2\)
\(a[1] = -e^{-2\pi f_c}\)

Lo que resulta en la siguiente respuesta en frecuencia

\[ H[k] = \frac{1+e^{-2\pi f_c}}{2} \frac{(e^{j \frac{2\pi}{N} k} - 1)}{(e^{j \frac{2\pi}{N} k} - e^{-2\pi f_c})} \]

Es decir un cero en \(1\), un polo en \(e^{-2\pi f_c}\) y ganancia \(\frac{1+e^{-2\pi f_c}}{2}\)

9.2.3. Aplicar un filtro a una señal con scipy¶

Para filtrar una señal unidimensional con un filtro IIR (sin variar la fase de la señal) podemos utilizar la función

    scipy.signal.filtfilt(b, # Coeficientes del numerador
                          a, # Coeficientes del denominador
                          x, # Señal a filtrar
                          ...
                         )

Los siguientes ejemplos muestran un señal de tipo pulso rectangular filtrada con sistemas IIR de primer orden pasa bajo y pasa-alto diseñados con las recetas mostradas anteriormente

n = np.arange(0, 500)
x = 0.5 + 0.5*scipy.signal.square((n)/(2.*np.pi*5), duty=0.3)

def iir_low_pass(signal, fc):
    gamma = np.exp(-2*np.pi*(fc))
    b, a = [(1-gamma), 0], [1, -gamma] 
    return scipy.signal.filtfilt(b, a, signal)

y = {}
for fc in [0.05, 0.02, 0.01]:
    y[fc] = iir_low_pass(x, fc)

px = hv.Curve((n, x))
py = []
for fc, y_ in y.items():
    py.append(hv.Curve((n, y_), label=f'fc={fc}'))

hv.Layout([px, hv.Overlay(py)]).cols(1).opts(hv.opts.Curve(height=200))

def iir_high_pass(signal, fc):
    gamma = np.exp(-2*np.pi*(fc))
    b, a = [(1+gamma)/2, -(1+gamma)/2], [1, -gamma]
    return scipy.signal.filtfilt(b, a, signal)

y = {}
for fc in [0.01, 0.02, 0.05]:
    y[fc] = iir_high_pass(x, fc)

px = hv.Curve((n, x))
py = []
for fc, y_ in y.items():
    py.append(hv.Curve((n, y_), label=f'fc={fc}'))

hv.Layout([px, hv.Overlay(py)]).cols(1).opts(hv.opts.Curve(height=200))

Nota

El filtro pasa-bajos suaviza los cambios de los pulsos rectangulares. El filtro pasa-altos elimina las zonas constantes y resalta los cambios de la señal.

9.3. Diseño de filtros IIR de segundo orden¶

Los filtros IIR de segundo orden o biquad tienen dos polos y dos ceros.

Su respuesta en frecuencia es

\[ H[k] = \frac{b[0] + b[1] W_N^k + b[2] W_N^{2k}}{1 + a[1] W_N^k + a[2] W_N^{2k}} = K \frac{(W_N^{-k} - \beta_1) (W_N^{-k} - \beta_2)}{(W_N^{-k} - \alpha_1)(W_N^{-k} - \alpha_2)}, \]

donde \(W_N = e^{-j \frac{2 \pi}{N}}\) y la relación entreo coeficientes y polos/ceros es:

\[ b[0] = K, \quad b[1] = -K (\beta_1 + \beta_2), \quad b[2]= K \beta_1\beta_2 \]

\[ a[1] = - (\alpha_1 + \alpha_2), \quad a[2]=\alpha_1 \alpha_2 \]

Con arquitecturas de segundo orden se pueden crear filtros pasabanda y rechaza banda

9.4. Diseño de filtros IIR de orden mayor¶

Para crear los coeficientes de filtro IIR de orden mayor podemos usar la función

scipy.signal.iirfilter(N, # Orden del filtro
                       Wn,  # Frecuencias de corte (normalizadas en [0,1])
                       fs, # Frecuencia de muestreo
                       btype='bandpass', # Tipo de filtro: 'bandpass', 'lowpass', 'highpass', 'bandstop'
                       ftype='butter', # Familia del filtro: 'butter', 'ellip', 'cheby1', 'cheby2', 'bessel'
                       output='ba', # Retornar coeficientes
                       ...
                       )

El filtro Butterworth es óptimo en el sentido de tener la banda de paso lo más plana posible.

Otros filtros se diseñaron con otras consideraciones.

Los filtros IIR digitales están basados en los filtros IIR analógicos.

Observe como al aumentar el orden el filtro pasabajo IIR comienza a cortar de forma más abrupta

Hk = {}
for order in [1, 2, 5, 20]:
    b, a = scipy.signal.iirfilter(N=order, Wn=0.2, fs=1,
                                  ftype='butter', btype='lowpass', output='ba')
    freq, response = scipy.signal.freqz(b, a, fs=1)
    Hk[order] = np.abs(response)

p = []
for order, response in Hk.items():
    p.append(hv.Curve((freq, response), 'Frecuencia', 'Respuesta', label=f'orden={order}'))
hv.Overlay(p)

9.5. Comparación de la respuesta en frecuencia de filtros FIR e IIR del orden equivalente¶

Comparemos la respuesta en frecuencia de un filtro IIR y otro FIR ambos pasa-bajo con 20 coeficientes

Fs = 1
fc = 0.25
h = scipy.signal.firwin(numtaps=20, cutoff=fc, pass_zero=True, window='hann', fs=Fs)
b, a = scipy.signal.iirfilter(N=9, Wn=fc, fs=Fs, ftype='butter', btype='lowpass')
display(len(h), len(b)+len(a))

freq_fir, response_fir = scipy.signal.freqz(h, 1, fs=Fs)
freq_iir, response_iir = scipy.signal.freqz(b, a, fs=Fs)

p1 = hv.Curve((freq_fir, np.abs(response_fir)), 'Frecuencia', 'Respuesta', label='FIR')
p2 = hv.Curve((freq_iir, np.abs(response_iir)), 'Frecuencia', 'Respuesta', label='IIR')
hv.Overlay([p1, p2])*hv.VLine(fc).opts(color='k', alpha=0.5)

La linea negra marca la ubicación de la frecuencia de corte

Nota

El filtro IIR es mucho más abrupto, es decir filtra mejor, que el filtro FIR equivalente

Una desventaja del filtro IIR es que por definición introduce una desfase no constante en la señal de salida

freq_fir, delay_fir = scipy.signal.group_delay(system=(h, 1), fs=Fs)
freq_iir, delay_iir = scipy.signal.group_delay(system=(b, a), fs=Fs)

/home/phuijse/.conda/envs/info320/lib/python3.9/site-packages/scipy/signal/filter_design.py:688: UserWarning: The group delay is singular at frequencies [3.105, 3.111, 3.117, 3.123, 3.129, 3.135], setting to 0
  warnings.warn(

p1 = hv.Curve((freq_fir, delay_fir), 'Frecuencia', 'Desfase', label='FIR')
p2 = hv.Curve((freq_iir, delay_iir), 'Frecuencia', 'Desfase', label='IIR')
hv.Overlay([p1, p2])*hv.VLine(fc).opts(color='k', alpha=0.5)

¿Cómo se ve una señal filtrada donde se preserva la fase versus una donde no se preserva la fase?

Consideremos la señal rectangular anterior y apliquemos un filtro pasa-bajo IIR de orden 1

Esta vez compararemos el filtro con la función scipy.signal.lfilter y la función scipy.signal.filtfilt. La primera no preserva la fase mientras que la segunda si lo hace

Fs = 1
fc = 0.01
n = np.arange(0, 500)
x = 0.5 + 0.5*scipy.signal.square((n)/(2.*np.pi*5), duty=0.3)

b, a = scipy.signal.iirfilter(N=1, Wn=fc, fs=Fs, ftype='butter', btype='lowpass')
# No se preserva la fase
y_lfilter = scipy.signal.lfilter(b, a, x)
# Se preserva la fase
y_filtfilt = scipy.signal.filtfilt(b, a, x)

px = hv.Curve((n, x), 'Tiempo', 'Entrada')
py = []
py.append(hv.Curve((n, y_filtfilt), 'Tiempo', 'Salida', label=f'Fase constante'))
py.append(hv.Curve((n, y_lfilter), 'Tiempo', 'Salida', label=f'Fase no constante'))

hv.Layout([px, hv.Overlay(py)]).cols(1).opts(hv.opts.Curve(height=200))

Nota

En el caso donde no se preserva la fase podemos notar que la señal de salida está desplazada con respecto a la original. Además los cambios tienen una transición asimétrica

La función scipy.signal.filtfilt “arregla” el problema del desfase filtrando la señal dos veces. La primera vez se filtra hacia adelante en el tiempo y la segunda vez hacia atrás. Por ende no se puede aplicar en un escenario de tipo streaming donde los datos van llegando de forma causal.

En una aplicación causal donde se necesite preservar la fase debemos usar un filtro FIR.

9.6. Apéndice: Efectos de audio con filtros IIR¶

El siguiente ejemplo muestra como implementar el conocido filtro Wah-wah usando un sistema IIR

Este es un filtro pasabanda modulado con ancho de pasada fijo \(f_b\) [Hz] y una frecuencia central variable \(f_c\) [Hz], donde La frecuencia central se modula con una onda lenta

Se modela como el siguiente sistema IIR

\[ H[k] = \frac{(1+c)W_N^{2k} -(1+c) }{W_N^{2k} + d(1-c)W_N^k -c} \]

donde

\[ d=-\cos(2\pi f_c/f_s) \]

y

\[ c = \frac{\tan(\pi f_b/f_s) -1}{\tan(2\pi f_b /f_s)+1} \]

Veamos como modifica este filtro una señal de audio

import librosa
data, fs = librosa.load("../../data/DPSAU.ogg")
Audio(data, rate=fs)

data_wah = []
zi = np.zeros(shape=(2,))
# Parámetros fijos del filtro
fb, Nw = 200, 5 
c  = (np.tan(np.pi*fb/fs) - 1.)/(np.tan(2*np.pi*fb/fs) +1)
# Filtramos una ventana de la señal moviendo lentamente fc
for k in range(len(data)//Nw):
    # Cálculo de la frecuencia central
    fc = 500 + 2000*(np.cos(2.0*np.pi*k*30./fs) +1)/2
    d = -np.cos(2*np.pi*fc/fs)
    # Coeficientes del filtro
    b, a = [(1+c), 0, -(1+c)], [1, d*(1-c), -c]
    # Filtramos, usando el filtrado anterior como borde (zi)
    data2, zi = scipy.signal.lfilter(b, a, data[k*Nw:(k+1)*Nw], zi=zi)
    # Guardamos
    data_wah.append(data2)

Audio(np.hstack(data_wah), rate=int(fs))

Si quieres profundizar en el tema de los filtros IIR aplicados a efectos de audio recomiendo: https://www.ee.columbia.edu/~ronw/adst-spring2010/lectures/lecture2.pdf

INFO183 Análisis de sistemas lineales

9. Diseño de sistemas y filtros IIR¶

9.1. Definición de un sistema IIR¶

9.1.1. Respuesta en frecuencia del sistema IIR¶

9.1.2. Ejemplo de respuesta al impulso de un sistema IIR¶

9.1.3. Ejemplo de respuesta en frecuencia de un sistema IIR¶

9.2. Diseño de filtros IIR simples¶

9.2.1. Receta para un filtro pasa bajo IIR con frecuencia de corte \(f_c\)¶

9.2.2. Receta para un filtro pasa alto IIR con frecuencia de corte \(f_c\)¶

9.2.3. Aplicar un filtro a una señal con scipy¶

9.3. Diseño de filtros IIR de segundo orden¶

9.4. Diseño de filtros IIR de orden mayor¶

9.5. Comparación de la respuesta en frecuencia de filtros FIR e IIR del orden equivalente¶

9.6. Apéndice: Efectos de audio con filtros IIR¶